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Y ■ Abstract 

pi i' We derive the long-wavelength effective action for the collective modes in sys- 

^ , terns of fermions interacting via a short-range s-wavc attraction, featuring un- 

equal chemical potentials for the two fermionic species (asymmetric systems). 
"j3 , As a consequence of the attractive interaction, fermions form a condensate that 

spontaneously breaks the t/(l) symmetry associated with total number con- 
servation. Therefore at sufficiently small temperatures and asymmetries, the 
'^ ' system is a superfiuid. We reproduce previous results for the stability condi- 

^ , tions of the system as a function of the four-fermion coupling and asymmetry. 

We obtain these results analyzing the coefficients of the low energy effective La- 
grangian of the modes describing fluctuations in the magnitude (Higgs mode) 
and in the phase (Nambu-Goldstone, or Anderson-Bogoliubov, mode) of the 
kJ , difermion condensate. We find that for certain values of parameters, the mass 

^sg ' of the Higgs mode decreases with increasing mismatch between the chemical 

OO . potentials of the two populations, if we keep the scattering length and the gap 

C^^ ' parameter constant. Furthermore, we find that the energy cost for creating a po- 

^^ \ sition dependent fluctuation of the condensate is constant in the gapped region 

Tlj" ■ and increases in the gapless region. These two features may lead to experimen- 

^D \ tally detectable effects. As an example, we argue that if the superfluid is put in 

rotation, the square of the radius of the outer core of a vortex should sharply 
increase on increasing the asymmetry, when we pass through the relevant region 
in the gapless superfluid phase. Finally, by gauging the global t/(l) symmetry, 
we relate the coefficients of the effective Lagrangian of the Nambu-Goldstone 
;J] ■ mode with the screening masses of the gauge field. 
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1. Introduction 

Experiments with trapped cold atomic gases have driven a renewed in- 
terest in fermionic pairing [ij, y]- In particular, much effort has been de- 
voted to understanding the superfluid phases of imbalanced fermionic gases, 
featuring unequal number of particles of the distinct fermionic species that 

pair d i, i, S, 0, i, i, [13, El, E E, [3 [II [II El M, M, M MM,M,^ I ■ 



The system consists of fermions of two different species, -01 and ip2, which 
correspond to two hyperfine states of a fermionic atom like ^Li. These fermions 
have opposite spin and the interaction between them can be tuned by employing 
a Feshbach resonance [25| . The strength of the interaction is given in terms of 
the s-wave scattering length between the two species. 

For zero imbalance, the system properties are qualitatively well understood 
using mean field theory [26|. In weak coupling the system lives in a weakly 
coupled BCS state and crosses over to a strongly coupled BEC state through 
the resonance region. While the extreme BCS and BEC regimes are also in 
good quantitative control in mean field theory, close to resonance (the unitarity 
region) a quantitative understanding of the phases comes mainly from Monte- 
Carlo calculations 11 1. (For other approaches see [27|, |29|, |28|.) This is because 



close to resonance the scattering length is much larger than the inter-particle 
distance and there is no small parameter in the Lagrangian to expand in. There- 
fore fluctuations may change the mean field results substantially. 

In standard BCS superfiuids the chemical potentials of the two fermionic 
species are equal. An imbalance in the number of V'l and ^'2 is implemented by 
taking the chemical potentials for the two species, /ii and fi2 respectively, to be 
different. (We will name our species in a way that /ii > fi2.) If the chemical 
potential difference, 2dfi ~ fJ.i — fJ-2 is much smaller than the magnitude of the 
gap parameter |A|, the splitting cannot disrupt BCS superfluidity because the 
superfluid state with equal number densities is energetically favored in compar- 
ison with a normal state with a fermionic imbalance. On the other hand, as 
pointed out in [3|, in the weak coupling regime, BCS superfluidity cannot persist 
for large values of S^. Indeed, there exists an upper limit for S^ (the so-called 
Chandrasekhar-Clogston limit), beyond which the homogeneous superfluid state 
is no longer energetically favored over the normal phase. 

For imbalanced systems, a qualitatively complete picture of the phase dia- 
gram has not been established yet. Proposed possibilities are phase-separation 
[7[ , breached pair superfluidity [J, |8|,j9|, [lO[ , deformed Fermi sea pairing |6| and 
non-homogeneous or LOFF pairing [5J. (See [30| and [31| for reviews.) 

The phase diagram of the system at T = as a function of the scattering 
length and the chemic al p otential difference has been explored in the mean 
field approximation in [23|, [l3|, [iGJ, [S^- The authors find that on the BCS side 



of the resonance there are no stable homogeneous superfluid phases that have 
gapless Fermi surfaces. On the BEC side of the resonance, there are stable 
gapless superfluid phases, which can exhibit a net polarization. At resonance, 
mean field theory suggests a first order phase transition from the superfluid to 
the normal phase as Sfi is increased, without any intervening gapless superfluid 



phase. Consequences of the phase diagram for experiments with trapped atoms 
were explored in [iGJ, llSJ . At resonance if we fill different number of -01 and '!/'2 in 
the harmonic trap, because the gapped phase can not feature a net polarization, 
the system phase separates with an unpolarized superfluid in the central region 
of the trap and a polarized normal fluid at the exterior. 

For non-zero imbalance close to the resonance, fluctuations may change the 
mean field results qualitatively. This has to be contrasted with the zero imbal- 
ance case, where fluctuations lead only to a quantitative change of the mean 
field results. Indeed for non-zero imbalance many features of the phase diagram 
are not caught by the mean field approximation. The authors of |15| g o beyond 
mean field theory by using results from Monte-Carlo simulations [2^ and pro- 
pose a phase diagram which features a splitting point near resonance at non-zero 
(5/x, where the homogeneous superfluid, a LOFF like inhomogeneous phase, and 
the gapless superfluid phase coexist. They also find stable gapless fermionic 
modes with one and two Fermi surfaces, on the BCS side of the resonance. A 
detailed treatment of fluctuations around the resonance using an expansion in 



e = D — A space dimensions at T = [33|, l29| supports this picture. A different 



approach consists in generalizing the Fermi gas to a model with 27V hyperfine 
states, performing a systematic 1/A'' loop expansion around the BEC-BCS so- 



lution [3J, |27| . The phase diagram at unitarity has also been explored using 
a Superfiuid Local Density Approximation (SLDA) [3^, [3g|. With this method 
one finds that on increasing Sfx from zero at unitarity, there is an intervening 
window of values for which the LOFF phase is favored over the homogeneous 
superfiuid and the normal phases. 

In this paper, we study small fluctuations about the mean field value of 
the gap parameter for a system with mismatched Fermi surfaces. We consider 
fluctuations of A both in its phase and in its magnitude. Both of these involve a 
coherent change in the wavefunctions of fcrmions in many different momentum 
eigenstates, and are therefore collective modes of the system. In particular, long 
wavelength fluctuations in the Nambu-Goldstone field are associated with the 
hydrodynamic mode (or sound mode) in the paired system and can be related to 



dynamic phenomena like compressions in a trapped atomic gas [38[ . By looking 
at the stability of the energy with respect to these excitations, we can map 
out the parameter values for which BCS-like pairing is favoured. We also use 
the expansion of the free energy in the magnitude of A to explore the typical 
length scale of inhmogeneities in the condensate in non-uniform configurations 
like vortices. 

2. Methods and materials 

In our quest to understand how fluctuations in the condensate about the 
mean field value affect the phase diagram of cold atomic gases with unequal 
number of ipi and 'ip2 fermions, we study the effective Lagrangian density de- 
scribing these fiuctuations. We do this by integrating out the fermions from the 
system and writing the effective action as a series in powers of the fiuctuations 



and their derivatives [37|, |28[ . We explicitly calculate the terms up to second 



order in the fluctuations and their derivatives. We expect that our mean field 
calculation of these coefficients are under better control away from unitarity [38| . 

For zero imbalance, the collective modes associated with fluctuations in the 
phase and the magnitude of the condensate were analyzed over the full BCS- 
BEC crossover in [39J . In the limit of long wavelengths (or small momenta) the 
theory is dominated by the Nambu-Goldstone mode associated with the phase 
fluctuations, travelling with the speed of sound given by c^ = {n / ra){d^ / dn) . 
The study by [39(] writes the effective Lagrangian to all order in derivatives, but 
only to the second order in fields. Very recently, in [40[ the effective Lagrangian 
describing interaction terms between the Nambu-Goldstone mode and the Higgs 
mode were obtained. By integra ting out the Higgs mode, the expression of the 
speed of sound first obtained in [4l| was reproduced in [40[ ■ 

In our study we restrict ourselves to only terms upto the second order in 
a derivative expansion. We reproduce the results of [39|, |40j and extend the 
analysis to non-zero imbalance. This is a physically interesting case because 
experiments have been performed for unequal number of ipi and il)2t and a 
change in the behavior of the collective modes can possibly give us information 
about novel phases that may arise in these experiments. In particular, we find 
that the Higgs mode mass shows an intersting behavior in the gapless BEG 
region as we discuss below. Because of this, we do not integrate out the Higgs 
mode as done by [40| , and keep it in the effective Lagrangian. 

Efforts to study the collective modes beyond the mean field approximation, 
by methods that may be under better control near unitarity, can be found 

[mil, Em. 
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The coefficients of the terms in the effective action tell us about the stability 
of the mean field solutions. The analysis of the stability of various phases in 
imbalanced Fermi gases has been studied previously in several different works. 
In [l6[ the authors looked at the phase diagram in detail, both in the narrow 
and the broad resonance limits. One important conclusion from their study is 
that it is important to check that the free energy is a local minimum rather 
than a local maximum, at the solution of the gap equation. In Ref. [14| it is 
shown that this criterion is equivalent to the requirement that the number sus- 
ceptibility is positive. In terms of the coefficients in the effective Lagrangian, 
it corresponds to the requirement that the mass-squared of the Higgs field be 
positive, ensuring stability with respect to homogeneous fluctuations. The au- 
thors of [l6| also derived the Ginzburg-Landau theory in the BEC regime for 
imbalanced Fermi gases, upto A^ in the fluctuations for the Higgs fleld about 
the normal phase (A = 0). The motivation for considering a Ginzburg-Landau 
expansion is that the gap is zero in the normal phase and expected to be small 
close to the gapless superfluid-normal phase boundary. This Ginzburg-Landau 
expansion can therefore be used to map the phase boundary between the two 
phases [l6[. Our expression for the quadratic coefficient in a Ginzburg-Landau 
expansion (shown in [Appendix G\ can not be directly compared to the expres- 
sion in [16| since this specific expression was given only in the narrow resonance 
approximation, while we work in the mean field approximation. However, by 
considering the stability of the Higgs field we conclude that there are locally 



stable gapless phases in the BEC regime, which go to the normal phase as we 
increase Sfi. This conclusion matches the conclusion by [16[ . We also go further 
by looking at the lowest non-trivial terms in the derivatives of the Higgs field. 

Several groups have analyzed stability with respect to space-time depen- 
dent (inhomogeneous) fluctuations in the condensate. This can give additional 
information to the study of local and global instability of homogeneous conden- 
sates because a phase could be stable with respect to a homogeneous change 
in the order parameter, but could be unstable with respect to the formation of 
inhomogeneous condensates. 

More specifically, the instability of gapless states towards the growth of phase 
modulation of the condensate, the so called current instability, has been studied 
in [isl [2^ |l9[ . The instability towards a growth of change in the magnitude 
(which will happen if the Higgs mass is imaginary), the so called Higgs insta- 
bility, was also studied by the authors of [19| . They showed that the absence of 
this instability is equivalent to the requirement that the number susceptibility 



matrix is positive definite J13l |22| . It was also found that the current instability 
is much less stringent than the Higgs instability. (For the manifestation of the 
current instability in the context of pairing in quark matter, see |44l . l45l . Uql . The 
Higgs instability in the quark matter context has been studied in [47|, |48|, |49| ) . 

By looking at the constraints on the positivity of the coefficients of the 
effective action, we reproduce the above mentioned results for stability. In 
addition, we consider the implication of the requirement that the energy cost of 
creating a position dependent fluctuation in the magnitude of the condensate 
(Higgs elasticity) be positive. This criterion has not been analyzed before in the 
literature, but we find that it gives a weaker condition than current stability. 

The main new results in the present paper are related to a study of the 
Higgs mass and Higgs elasticity as a function of the coupling and the chemical 
potentials. We find that the Higgs mass is small in the gapless phase in the 
BEC regime. We also find that in the BEC regime, the Higgs elasticity is 
constant in the gapped phase and increases in the gapless phase. This has 
important consequences for any non-homogeneous configuration created in a 
system tuned to sit in this region. It implies that a cost of creating a gradient 
in the condensate value is large, and hence the condensate should vary slowly in 
any such configuration. An inhomogeneous configuration has been considered 
by [50| , who however evaluaed the elasticity of the condensate field at unitarity 
for vanishing values of the gap, while we consider fluctuations about the mean 
field solution. Our analysis is also an improvement over the analysis of j5ll | 
where the Higgs elasticity is not computed microscopically. 

This paper is organized as follows. In Section I3TT1 we present our model and 
review some basic equations of the mean field analysis. In Section F3.2[ we study 
the fluctuations of the difermion condensate, and derive the general expression 
for the effective action for the fluctuation fields, valid up to second order in the 
fluctuations. We consider both fluctuations in the magnitude, and in the phase 
of the difermion condensate. The reader not interested in the calculational 
details may skip over to Section l3.3i where we present the low energy effective 
theory for these modes. We show the expressions of the coefficients that appear 



in the effective Lagrangian for arbitrary values of the temperature. From the 
sign of these coefficients we obtain stabihty criteria that we analyze in detail 
in the case of vanishing temperature. From this analysis we reproduce the 
conclusion that there exist stable gaplcss phases in the BEC (strong coupling) 
regime at non-zero asymmetry. The central results of the paper are discussed 
in Section 14.11 where we evaluate the mass of the Higgs mode in the strong 
coupling regime and find that for certain values of the parameters in the gapless 
region, this mode is light. This implies that the outer core of a vortex in this 
region will be wider than in the gapped superfiuid phase. Furthermore, we find 
that the elasticity of the Higgs mode sharply increases in the gapless region on 
increasing the asymmetry. This also means that the radius of vortices will be 
large in this region. This effect could be experimentally detectable in cold atoms 
experiments where the mismatch between the two species can be tuned. 

In [Appendix A| wc show the equivalence between the coefficients in the 
effective action of the Nambu-Goldstonc mode and the screening masses that 
are obtained by gauging the U{1) symmetry. Formally, the equivalence may 
seem apparent from gauge invariance, but the explicit demonstration of the 
same is non trivial, and therefore we include the derivation in [Appendix A 



In [Appendix B[ we report some details of the calculation of the coefficients 
appearing in the effective Lagrangian. 

3. Calculation 

3.1. Model and ansatz 

Wc consider a non-relativistic system consisting of two species of fermions 
f/^i and '02 of equal mass m but different chemical potentials /ii = fi + S^ and 
/X2 = /x — (5^, with /J, being the average of the two chemical potentials and 2S^ 
the difference between them. Defining the field tp = {ipi ip2)'^ , the Lagrangian 
density describing free fermions can be written as, 

Cf^ij^{idt~Eip)+fi + 6fia'')ij, (1) 

where E{p) — p^/(2to), with p the momentum operator V/i. The energy of a 
free fermion relative to the average chemical potential is conventionally indicated 
by 6(p) = E{p) — ^. We assume that the Feshbach interaction between fermions 
of different species can be modeled by a point like four Fermi interaction, and 
the corresponding term in the Lagrangian can be written as 

^i = ^i^Ux)^lix)M^)M^) , (2) 

with a,/3 G {1,2} and where A > for attractive interaction, the case we are 
interested in. 

The effect of the attractive interaction between fermions is to produce a 

difermion condensate 

A(x) 

{'Ipaix)'ipi3{x)) ^ -^EafS , (3) 



where e is the two dimensional antisymmetric tensor e = ia'^ . 

In the mean-field approximation the Lagrangian density can be written as 

|A(x)p 



£ = *^ 



idt - C(P) + Sn<y^ 
A*(a;)£ 



-A(a;)e 



* 



A 



(4) 



where ^ stands for the four component Nambu-Gorkov spinor, 



* 



1 

71 






(5) 



The fluctuations of the condensate will be treated in the next Section. Here we 
only discuss the homogeneous phase, with A(a;) = A = const. In this case the 
excitation spectrum is described by the quasiparticle dispersion laws 



e+ = +(5Ai+Ve(p) +A 



e_ 



5/i + Ve(p) +A' 



(6) 



The knowledge of the dispersion laws of the system allows one to evaluate the 
grand-potential, which is given at T = by the expression. 



_ A^ 1 






|e+| + |e_|-2e(p) 



(7) 



The integral in this expression is ultraviolet divergent and can be regularized in 
the usual way j52l | , by writing A in terms of the scattering length a according to 



m 
A'Ka 



1 
A 



cPp 1 
(27r)3^ 



For later convenience we introduce the dimcnsionlcss coupling constant 



9 



1 
kpa 



(8) 



(9) 



where kp is the Fermi momentum of the system which is defined in terms of 
the average number density n of the two species by the relation n = kp/{3Tr'^). 
The weak coupling regime, where the BCS approximation holds, corresponds to 
(7 — >■ — oo. This approximation is generally very good for superconductivity in 
metals. On the other hand, in cold atoms the strength of the interaction can 
be varied in the vicinity of a Feshbach resonance, where the scattering length 
strongly depends on the applied magnetic field. Therefore both the weak and 
strong coupling regimes can be reached in this case. 

Knowing the free energy of the system, one can evaluate the gap parameter 
A by solving the equation 



Let us note explicitly that we do not write equations for /ii and fi2- We do 
not work at fixed particle number densities ni and 712 and therefore we do not 
impose the equations: 

dn on 

t: — = -"■! 7^ — == -"■2 , (11) 

which would be needed in the analysis if ni and ^2 were held fixed [l3|. Instead, 
the values of n's for given A and /i's can be determined by the relations, Eq. pTj) . 
Note also that the effect of the condensate is to spontaneously break the 
global C/(l) symmetry corresponding to the conservation of the total fermion 
number, ni +712. Therefore there will be a Nambu-Goldstone mode associated 
with the spontaneous breaking of this symmetry, and the system will conse- 
quently be a superfluid. Clearly if one gauges this symmetry, the spontaneous 
breaking of the local symmetry leads to the appearance of a mass term for the 
gauge boson (Meissner mass), and the system becomes a superconductor. In 
the following analysis we will assume that the U{1) symmetry is global, i.e. 
fermions are not charged, and therefore we will study the dynamics of the asso- 
ciated Nambu-Goldstone boson. In [Appendix A we will consider the relations 



between the parameters appearing in the Lagrangian describing the Nambu- 
Goldstone bosons, and the screening masses of the gauge field. 

3.2. Fluctuations 

In order to include fluctuations of the condensate, we introduce the field 
r]{x) that represents the deviation of the condensate from its mean field value. 
In the presence of fluctuations, A (a;) in Eq. ((4|) is given by 

A(x) = A + 77(x) , (12) 

where it is assumed that the fluctuation is much smaller than A. In this pa- 
per, we will consider only homogeneous condensates, meaning that A on the 
right hand side of Eq. ()12p is independent of x. In principle one might con- 
sider the case where the underlying condensate is x dependent, like in the non- 
homogeneous LOFF phase. However, we will postpone the study of such a case 
to future work. In order to simplify the analysis, but without lack of generality, 
we choose the phase of the fermion flelds so that the mean field condensate, 
A, is real. The field ry, on the other hand will have both real and imaginary 
components. 

For a given temperature T, the partition function is given by. 



Z= I?77*2?772?*^2?1'e-'5[**'*''''''*1 , (13) 

where S is the Wick rotated action, 
S[i'\'i',ij,rf]^Jd^x[j\A + rj{x)\' 



|2 



-d,4-^{p)+dna^ -(A + ?7(x))e 



iA + 7f{x))e -d^. + Up) - Sfia' 



(14) 



and we use the imaginary time formalism where x^ is the imaginary time it, 
and runs from -1/{2T) to 1/(2T). 

To find the effective action for the rj field, we integrate out the fermionic 
field, which can be done because the action is quadratic in ^P. This gives. 



P?7*P?/e-'5["'"'l , (15) 



with 



- {^Trlog 



dx'i - C(p) + <5Ai 


-(A + ,7(x)) ^ 

-9,4 + e(p) + s^i ^ 


+ {6fi- 


(16) 



where Tr symbolizes the trace over Nambu-Gorkov indices and over a complete 
set of functions over space-time. The factor of 1/2 before the Tr takes care of 
the fictitious doubling of degrees of freedom that arose when we introduced the 
Nambu-Gorkov spinor. 

At a formal level, Eq. P^ gives the desired effective action for the fluctua- 
tions. However, it is not possible to compute the Tr analytically for arbitrary 
functions ri{x) and hence we expand the logarithm in increasing powers of rj 
(and 77*), 

Trlog(d + V)^ Trlog(O) + Tr(^ — (-Q-^f )") , (17) 

n—1 

where we have defined 

A ( A{p) -A\ ., _ I ( A(j>) A \ 

^ ^ ( -.'(.) 'T ) ' ^^^) 

and where 

Aip) = ipi - e(p) + 5^1 , Aip) = ipi + e(p) + Sn , (20) 

with p = (— 9,^, V/i) the (Euclidean) four momentum operator. The quantity 
appearing in the denominator of Eq. (jTSl) is given by 

Dip) = Aip)Aip) - A2 = {tpi + 5^i + e(p)) {ip^ + 5^l- e(p)) , (21) 



\2 : \2 



where we have also defined e(p) = v C(p) + A^ 

We thus obtain the effective action as a series expansion 

5[77,77*]=5("'+5(i)+5(2)+... , (22) 



with iS^'-* proportional to the zth power of rj (and 77*). We shall now analyze the 
various terms in this expansion individually. 

The zeroth order contribution to the action 5'°-' is proportional to the free 
energy of the system in the absence of fluctuations, 

p 

where V is the spatial volume of the system and fl is the free energy at finite 
temperature. In Eq. (j23|) and below, wc will use a notation where the sum over 
p refers to a sum over (spatial) momentum eigenvalues, p. and a sum over pi 
which runs over the fermionic Matsubara frequencies aj„ ~ {2n + 1)'!tT, for n 
integer. Bosonic Matsubara frequencies, w„ = 2m:T, will be denoted with k^ . 
Extremizing the free energy with respect to A we find two stationary points 
corresponding to the trivial solution A = 0, and 



lv^^)+i'^^^-'^) 



(24) 



In the following wc will assume that A is non-zero and use Eq. (|24p to simplify 
various expressions. 

We now turn to the term of the action in Eq. ()22p that is linear in 77, i.e. 
S^^' . This term is given by. 






where fj^k) and T]*{k) are the Fourier transforms oirj{x) and Tf{x) and are given 

by 



fl{k) = d'^xri{x)e'^'' 
rj*{k) = /d^Tr;*(x)e*'=-^ . (26) 



Employing the gap equation (Eq. ([M]) '). one obtains that 5^^-' = 0. This is 
clearly a consequence of the fact that we are considering a stationary point of 
the action. This result also holds if wc consider the solution A = 0. 

The lowest order non-trivial term in the expansion of the action is the one 



10 



quadratic in 77 

5(2) =jld^x {r?(x),f (x)} + iTr{(6-iy)2 + (Sf, ^ -S^,)} 

1 T v-r~ 



XV 

k 



Y.[fi{-kW{k)] 






2A{p)A{p + fc) 



fl{-k)yf{k) + {6fi -^ -SfJ.)\ , 



D{p)D{p + k) 

where the sum over k means integration over the three-momentum k and sum 
over bosonic Matsubara frequencies. 

Using the the gap equation we can simphfy the expression above as follows: 

5(2) = -|^,)(-fc>f*(fc){/2(fc) + 2/i(fc)}-^^,K-fc)'r(fc)/3(fc) 

k k 

-| E(^~*(-^)^~* W + fii-k)vik))hik) , (28) 

k 

where we have defined, 

^^(^) = 'TV^ Dip)Dip + k) +^'^'^-'^^ 

^,,, 1 T ^ {A{p + fc) - A{p)){A{p + fc) - A{p)) , ,^ 

^2(^) = iyE D{p)D{p + k) + is^.^-5^^) 

r ,,. -lT^ A{p)A{p + k)~A{p + k)A{p) , ^^ 

'^(^) = -v\ D{p)D{p + k) + (^^-^-^M). (29) 

Here Ii{k), l2{k) and Isik) are even in k; /i(fc) and /2(fc) are even in the 
time component k^ as well, while Isik) is odd in k^. Therefore we have that 
/i(— fe) = Ii{k), l2{—k) = l2{k), and /3(— fc) = —I^ik) with 4-d momentum fc. 
Note that the ultraviolet divergent contributions cancel exactly: /i, I2 and /s 
are all ultraviolet finite. 

In order to clarify the expression that we have obtained, it is convenient to 
separate r/ into its real and imaginary parts, 

Ti{x)^^{\{x)+id{x)). (30) 

Thus, the action 5'^' in Eq. (P5| can be written in terms of the A and 9 fields, 
as 



fe 



(31) 
11 



The evaluation of the functions Ii{k), l2{k) and laik) for arbitrary values of 
k is quite involved |28| . However, if we are interested in the long wavelength 
fluctuations of the condensate, we can expand the integrals in a power series in 
k and obtain the low energy effective action of the system. 

Note that for certain values of (5/z and A, the system may feature gapless 
fermionic modes that also contribute to the low energy dynamics of the sys- 
tem (H. 

In the following Section we will study the low energy effective Lagrangian of 
the system discarding the possible contribution of gapless fermions. This will 
allow to elucidate the role of the fields A and 6. 

3.3. Low energy effective Lagrangian 

The physical meaning of the real and complex components of the field ry, 
namely A and 6, is easy to understand in small fluctuation and long wavelength 
limit. We will show that in this limit A corresponds to the Higgs field and 9 to 
the Nambu-Goldstonc mode. 

Moreover in the limit of small k it is possible to expand /i(fc), /2(fc) and 
Iz{k) in a power series in k and to evaluate analytically or numerically each 
term of the expansion. 

Upon making this expansion, we obtain to second order in fc, 

/2(fc) = Akl-^\^ + 0{k^) 

/2(fc)+4/i(fc) = -C + Dkl-^]<?+0[k^) 

Ls{k) = ~koF + 0{k^), (32) 

where the expressions of the coefficients A, B, C, D, E and F are reported 
in [Appendix B[ As a check of our results we notice that taking 6ii ~ in. 
the expressions above, we reproduce the coefficient of the effective Lagrangian 
obtained in [40|. In particular we notice that for vanishing mismatch one has 
that A = AC, which matches with the result of [40|. However, for Sji =/= such 
a relation does not hold. 

To understand the physical meaning of the various coefficients in the effective 
action, let us first consider the case where the phase of the condensate, but not 
its magnitude, fiuctuates. That is, 

A -> Ae*'*^^) . (33) 

The field (j){x) represents the Nambu-Goldstone mode associated with the spon- 
taneous symmetry breaking of the total fermion number, jii + Ji2. Since there is 
no term that explicitly breaks this symmetry, the mass of this Nambu-Goldstone 
boson is exactly zero. 

Then the meaning of the coefficients A and B becomes clear if we notice that 
to linear order in 0, Eq. ((33|) corresponds to X{x) = and 6{x) = \/2A(j){x). 
The low energy Lagrangian density for is therefore, 

C^ = A^[A{dtcl>{x)r - |(a,(/.(x))2] . (34) 
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Therefore A is the coefficient appearing in the kinetic energy density of the 
Nambu-Goldstone mode, and B is related to the spatial variation of the Nambu- 
Goldstone mode. A negative value of i? or of A tells us that the mean field 
solution of the system is unstable to the growth of phase fluctuations of the 
condensate. If both B and A are positive, the system is stable. In our analysis 
we never find negative values of A, but we find regions of the parameter space 
where B is negative. 

The speed of sound, or equivalently the speed of the Nambu-Goldstone mode, 
is the same as the speed of cj) field, y^B/{3A), in the weak coupling BCS regime 
where integrating out the Higgs mode does not change the speed significantly. 
Therefore we reproduce the well known weak coupling result. It is easy to see 
that by integrating out the Higgs mode we reproduce the speed of the Nambu- 



Goldstone mode calculated by [38| and verified by [40|. Actually, from our 
expressions one can further extend their result to non-zero dfi. 

Now consider the case where only the magnitude of the condensate fluctu- 
ates, corresponding to the Higgs mode A ^- A + A(.t)/\/2. The low energy 
Lagrangian density for these fluctuations is, 

£, - -IcXixf + ^D{dtX{x)y- - |(a.A(x))2 , (35) 

that for positive values of the coefficients C, D and E is equivalent to the 
Lagrangian density of a massive bosonic field with mass squared (i.e. the square 
of the gap in the excitation spectrum) equal to C/D. If the various coefficients 
are not positive, then the system is unstable. We shall now analyze the three 
terms appearing in this Lagrangian. 

The CX^ term corresponds to the mass term and it can be interpreted as 
the change in the free energy, reported in Eq. (|23|) . caused by changing the 
magnitude of A. Since the mean field value of A is chosen so that the free 
energy is a local extremum, the sign of C tells us whether this extremum is 
a local maximum, for C < 0, or a local minimum, for C > 0. Therefore 
in the former case the system is unstable, in the latter it is stable or meta- 
stable, depending on whether the local minimum is also the absolute minimum 
of the system or not. It can be shown analytically (and we have also checked 
numerically) that the curvature of the potential around the stationary point is 
proportional to C. If one of the coefficient D or i;^ is negative, then the mean 
field value is unstable with respect to time- or space-dependent fluctuations of 
the magnitude of the condensate. We flnd that D is always positive, whereas E 
is negative in a certain region of parameter space. 

Negative values of B and E are both related to the growth of spatially non- 
uniform fluctuations of the condensate but may point to different possibilities 
for the true ground state of the system. A negative B may suggest that the 
condensate prefers to develop a non-zero phase modulation which carries a cur- 
rent, balanced by a counter-propagating current carried by gapless fermions. A 
non-zero E points to the formation of a spatial modulation in the magnitude of 
the condensate, which does not carry a current. 
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The coefficient F does not appear in the discussions of the Nambu-Goldstone 
and Higgs Lagrangians above. F mixes the A and 9 components. Such mixing 
between the components of a complex field has been discussed previously for 
Lagrangians featuring a global symmetry corresponding to phase rotations of 



the field 5J, |55| . In these cases the mixing term can be interpreted as a chemical 
potential for the conserved charge. Here, the two modes are not on an equal 
footing and the interpretation of this term may be more involved. We leave 
further discussion of the mixing between the Nambu-Goldstone and Higgs modes 
for future work. 

Note that in Eq. ([5^. the small momentum expansion of Ii{k), l2{k) and 
Isik) should be done with care. Namely, since in some cases these integrals 
are divergent, one cannot interchange the order of taking small k limit and p 
integration. 

3.4- Analysis of stability at T ~ 0. 

In the [Appendix B| we have reported the equations for the coefficients A, 
B, C, D and E for arbitrary values of the temperature. However, in the present 
paper we content ourselves with the analysis of the stability for the case of 
vanishing values of the temperature. 

Before considering the general case of arbitrary coupling, it is instructive to 
consider the limiting case of weak interaction. At weak coupling, the BCS hier- 
archy of scales, 6fj., A <C /i, holds. Therefore one can carry out the momentum 
integration analytically in a thin shell around the common Fermi surface, fi. 

Of particular interest is to study the phases which feature gapless fermionic 
excitations. These phases correspond to 6fj. > A and are known to be unstable 
in weak coupling. In this case, the coefficients appearing in the Lagrangian of 
the Nambu-Goldstone mode are given by 



1 TO(2m//)i/2 
8^ A2 



A = ^^^^^^i^H^.) 



1 (2m//)3/2 l-x 
while for the coefficients related to the Higgs mode we obtain 



D = 



1 TO(2m^)^/2 1-x^ 



3Sii 1 — X 



^ 1 (2m/.)3/2 l-x^ 

87r2 3mV x^l - x^) ' ^ ' 



where we have introduced x = \/ 5y? — A'^ j 8[i < 1. The mixing term is 

F = 0. (38) 
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The last equation shows that the Nambu-Goldstone and Higgs modes decouple 
in the weak coupling. Equations (p6)) and (p7)) show that both the Nambu- 
Goldstone and the Higgs fields develop instabilities in this regime, because the 
coefficients B, C and E are negative. The fact that B is negative indicates 
instability towards a phase with spontaneous generated currents [4J, |43, H, 47 1 
and a negative E towards a modulation of the magnitude of the condensate 48 , 
|49|. Negative C shows that this gaplcss phase does not correspond to a local 
minimum of the energy. However, in the weak coupling case it is known that 
well before the gapless phase develops, there is a first order phase transition 
to the normal phase or to a non-homogeneous superfluid phase. Indeed, for 
dn > A/v2 the energy of the local minimum corresponding to the non trivial 
solution of the gap equation is larger than the energy of the unpaired phase. This 
means that the Higgs and the Nambu-Goldstone modes that we are studying 
and that eventually become unstable at 6fi ^ A, correspond to fiuctuations 
around the meta-stable solution for Sfi > A/-\/2. 
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Figure 1: Regions in the (Sp,, p,) plane which feature zero, one or two spherical surfaces 
in momentum spaces with gapless excitations. These regions are marked with 0, I and II 
respectively and arc separated by solid lines. 



From this weak coupling analysis it is clear that in order to obtain a stable 
gapless state one should study the strong coupling regime realized for larger 
values of the coupling constant. To analyze the stability of the various super- 
conducting phases, we need to calculate the values of A for given values of A, 
fi and (5jU and then ascertain whether the coefficients A, B, C, D and E are 
positive. In particular, one of the questions we are interested in from such a 
study is whether there are regions of parameter space featuring stable phases 
having gapless excitations on one or two spherical surfaces in momentum space. 
One way to study this question without solving the gap equation is to eliminate 
the variable A by writing fi and J/i in units of A |22| . Therefore we define. 



^=A 






(39) 
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and the coefficients A, B, C, D, E and F are then functions of p, and Sp,, 
multiplied by appropriate powers of A and m to give the correct dimensions. 
We can then map out the region in the {5jl,jl) space where the integrals are 
negative, indicating instabilities. 

On the same {Sfl,p,) plane we can identify regions where the system has 
gapless excitations. Of the two dispersion laws reported in Eq. ^ the one 
indicated with e_ can become gapless in a certain range of parameters. This 
dispersion law is given by 

e_ (p) = -Sn + v/(pV(2m) - Ai)2 + A2 = A (-Sp + ^{p^-fiy + l) 

p (40) 



with p = 



^/2^^ ' 



and it can have zeros as a function of p (or p). 

In Fig. [1] we have divided the {Sp,, p) plane in three regions corresponding to 
the different number of gapless surfaces in momentum space and marked such 
regions with 0, I and II. The region marked with corresponds to (5/2 < 1 or 
5p > 1 with p < —y^Sp"^ — 1, where the dispersion law e_ has zero gapless 
modes. Region I corresponds to (5/2 > 1 and p e [—y^Spi'^ — 1, +y/Sp'^ — 1], 
where e_ is zero on one spherical surface in momentum space. Finally, the 
region II corresponds to Sp > 1 and p > +y^Sp'^ — 1 where e_ is zero for two 
distinct values of p, corresponding to two spherical surfaces in momentum space. 

In Ref. [24I an analysis of the stability of the various regions reported in this 
diagram has been done. In that paper the following requirements have been 
considered: 

i. The Meissner mass of two fictitious gauge bosons that couple to the 
fcrmions ■01 and -02 should be real and positive. 

ii. The 2x2 number susceptibility matrix associated with the two chemical 
potentials /ii and fi2 should be positive definite. 

iii. The free energy of the superconducting state should be lower than the 
free energy of the unpaired state, meaning that the pressure in the super- 
conducting phase has to be larger than the pressure in the normal phase. 

It turns out that the positivity of the Meissner mass leaves some region in 
the parameter space where the gapless state with two Fermi surfaces is stable. 
However, requiring the positivity of susceptibilities eliminates all the gapless 
states with two Fermi surfaces. Considering all the stability criteria above leaves 
only a narrow strip at /i < 0, where the gapless state with one Fermi surface is 
stable. 

We conduct a similar study by requiring that the coefficients A, B, C, D 
and E are positive. It turns out that A and D arc positive in all parameter 
space, while the other coefficients are negative in some regions. Since A and D 
turn out to be positive in the whole ((5/2, p) plane, the requirements of stability 
can be expressed throught the following criteria: 
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1. The coefficients B and E must be positive. This corresponds to have a real 
speed of sound for the Nambu-Goldstone mode and for the Higgs mode. 

2. The coefficient C must be positive. This corresponds to requiring that the 
superfiuid state is a local minimum of the free energy. 

3. The free energy of the superfiuid state should be lower than the free energy 
of the unpaired state, i.e. fls — fin < 0, where fig and i7„ refer to the free 
energies of the superfiuid and the normal phases, respectively. 

Notice that according to [1J| the stability criterion 2 is equivalent to criterion 
ii. above. The stability criteria 2 and 3 have been used to map out the phase 
diagram of imbalanced Fermi gases at both zero and non-zero temperatures in 
Refs. [ly, [ij, [ij, |22|, |23|. In these papers it is shown that the most stringent 



condition, for any values of /i and Sfi^ is that the free energy of the superfiuid 
state should be lower than the free energy of the unpaired state, corresponding 
to criterion 3 above. Here we want to remark that the requirement that the 
coefficients B, C and E arc positive, does still give some information about the 
system. Consider as an example gaplcss states that satisfy criteria 1 and 2, but 
fail 3. In this case the system is in a metastable gapless states that may be 
realized and studied in experiments. 

In Fig. [2] we report the results of our analysis concerning the stability criteria 
1, 2 and 3 above. On the left panel we report the results regarding the stability 
criteria 2 and 3. Criterion 3, corresponding to the requirement that ilg — r2„ > 0, 
excludes the shaded region directly above the blue dashed line. Criterion 2, 
corresponding to the requirement C > 0, excludes all the region directly above 
the dotted green line. A comparison with |22| shows that the requirement that 
C > is equivalent to the condition that the number susceptibility of the 
system should be positive. Therefore there is a sliver of parameter space where 
the superfiuid phase is meta-stablc and not absolutely stable. Our findings are 
in agreement with the results of Ref. [l6[, where it is found that deep in the 
BEG region a meta-stable gapless state exists. 

While criterion 1, corresponding to the requirement that B and E are pos- 
itive, is not as restrictive as criterion 2, it is still interesting because it tells us 
about the tendency of the system to turn into a non homogeneous phase. On 
the right panel of Fig. [2] wc report the results of the stability analysis concerning 
the coefficients B and E. The requirement B > excludes the shaded region 
of parameter space directly above the dotted green line and is equivalent to the 



requirement that the Mcissncr mass be real 22|. Requiring E > excludes 
the region directly above the dot-dashed red line. Therefore, the requirement 
i? > is more restrictive than the requirement E > 0. Hence wc find that 
the additional consideration of the position dependent fiuctuations in the Higgs 
field does not yield a more stringent criterion for stability than the requirement 
that there be no current instability. 

Notice that these criteria do not forbid the existence of states with two gap- 
less surfaces. For reference, the dshcd blue curve corresponding to the criterion 
3 is also reported in the right panel of Fig. [51 We now look at the implications 
of the variation of the expansion coefficients as a function of Sfi, fi and A, for 
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Figure 2: (color online) Left panel: Stable and unstable regions in the (5/i, /I) plane according 
to criteria 2 and 3. The requirement f!s — S7„ > 0, (criterion 3) excludes the shaded region 
directly above the dashed blue curve. We will refer this curve as "curve 3" below. The 
requirement C > (criterion 2), excludes the region directly above the green dotted curve. 
We will refer to this curve as "curve 2" below. Criterion 3 is more restrictive than criterion 2; 
there is a sliver of parameter space where the superfluid phase corresponds to a local minimum 
but not a global minimum of the free energy. Right panel: Regions in the (5/1, /I) plane which 
are stable or unstable according to the criterion 1 and 3. The requirement _B > excludes 
the shaded region directly above the dotted green line. The requirement i? > excludes the 
region directly above the dot-dashed red line. We see that the requirement B > is more 
restrictive than the requirement E > 0. In any case, these two requirements leave regions 
of parameter space showing two gapless surfaces, however this region is excluded once the 
criterion 3, corresponding to the dashed blue line, is considered. On the top of both figures 
the Chandrasekhar-Clogston limit 5/.t/A = l/\/2 ~ 0.707 is indicated, which corresponds to 
the critical value of the chemical potential splitting for the favorability of the superfluid phase 
in weak coupling. 0, I and II refer to the regions with zero, one and two gapless surfaces 
respectively, as in Fig. [T] 



the variation of the length scale of the modulation of the condensate in vortices. 

4. Results and discussion 

J^.l. Parameters of the Higgs Lagrangian and vortex radius 

The requirement that small fluctuations in the magnitude and the phase of 
the order parameter increase the free energy rather than decrease it, provides a 
strong constraint on the values that A, ^ and S^ can take in asymmetric cold 
atomic systems. The strongest constraint from these "local criteria" comes from 
the requirement that the value of A be a local minimum of the free energy rather 
than a local maximum (criterion 2). This condition excludes the possibility 
that there can be two spherical surfaces in momentum space featuring gapless 
quasiparticle excitations. A stronger constraint is provided by a global condition 
that the homogeneous superfluid phase has a lower free energy than the normal 
phase (criterion 3). 

From Fig. [21 one can notice that for Sfx > A, the curves associated with 
criterion 2 and criterion 3 run very close in the gapless region. Indeed these two 
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Figure 3: (color online) Curves of constant k in the Sfi, p, plane. Shown are four curves (dot- 
dashed lines, purple online) corresponding to k = —0.5, k = 0, k = 1 and k = 1.71. Negative 
scattering lengths, and therefore negative k, correspond to the BCS regime while positive k to 
the BEC regime. We look at the Higgs mass as a function of Sp, along the curve corresponding 
to K = 1.71 in Fig.U Point P, corresponding to Sp, ~ 1.59, has the largest value of Sp along 
the curve, for which Qs — f^n < 0. Between P and Q, corresponding to 1.59 < Sp < 1.66, the 
superfluid phase is metastable. To the right of the point Q, corresponding to Sp > 1.66, the 
superfluid phase is locally unstable, meaning C < 0. 



curves appear to converge asymptotically, for S^^ A. We recall that the Higgs 
mass is zero along curve 2 (dotted line (green online) in Fig. [2]). This suggests 
that the mass of the fluctuations in the magnitude of the condensate is very 
small along curve 3 (dashed line (blue online) in Fig. [2]) in the region /, and 
gets smaller as the two curves come closer. Since the presence of a light Higgs 
mode may be experimentally detectable, we have explicitly studied the mass of 
the Higgs field in the region where /i < 0, as a function of Sfi. This region in 
parameter space is accessible with positive values of the scattering length a, and 
lies on the BEC side of the resonance. 

To be concrete, we first solve the gap equation for various scattering lengths 
and see where we land in the parameter space. The result is shown in Fig. |31 
The four dot-dashed lines (purple online) show how /2 varies as a function of Sp, 
for four different values of the dimensionless variable, k = n / {2\/2mAa) [22|. 
Values of K ^ — 1 correspond to being deep in the BCS regime, while k ^ 1 
corresponds to being deep in the BEC regime. Since we are interested in the 
BEC regime we consider, for definiteness, the curve corresponding to k = 1.71. 
It intersects the curve corresponding to criterion 3 in P, at Sp, ~ 1.59, and the 
curve corresponding to criterion 2 in Q, at (5/2 ~ 1.66. 

Now consider the value of C/D, which is the mass squared of the Higgs 
fluctuation of the condensate, as we increase Sp along the curve labeled k, = 1.71 
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Figure 4: (color online) Mass squared of the Higgs mode, m,^ = C/ D, 
function of &fi along the curve k = 1.71 (see Fig. [Sj. In the gapped region m|j is a constant, 
and decreases when we enter the gapless regime. At Sp, ~ 1.59 the Higgs mass has the smallest 
value in the regime where the homogeneous superfluid is favored over the normal phase. For 
&fi > 1.66, corresponding to points on the right of Q in Fig. [S] the Higgs mass becomes 
imaginary and the superfluid phase is locally unstable. 



in Fig.[3l (See Fig. HI) As can be seen in Fig. SI for (5/2 = 0, the superfluid phase 
is favored over the normal phase and is also locally stable, meaning C > 0. As 
we increase 5p., as long as we are in the gapped phase, the free energy of the 
superfluid phase is independent of (5/2, (although ri„ decreases as we increase 
(5/.i) and hence the mass squared of the Higgs is positive and independent of 
(5/2 in this region. As we cross into the region featuring one gapless surface, C 
decreases as we move closer to the curve 2. When (5/2 ~ 1.59, corresponding 
to point P in Fig. [31 we have reached the largest value of Sjl for which the 
superfluid phase wins over the normal phase. This gives the smallest value of 
the Higgs mass in the region where it describes oscillations about the global 
minimum. We note that rri^^ drops by a factor of about 7.5 at (5/2 ~ 1.59 from 
its value at (5/2 = 0. 

Moving along into the metastable region between point P and Q, the Higgs 
mass square decreases and finally becomes negative when wc cross curve 2 at 
point Q in Fig. [3l corresponding to (5/2 ~ 1.66. 

Note that this calculation is done in a region where mean field methods are 
expected to be reliable. To illustrate this, we calculate the value of the inverse 
of the dimensionless expansion parameter, g = l/{kpa). Large and negative 
values of g correspond to being deep in the BCS regime while large and positive 
values of g correspond to the region deep in the BEC regime. At the point P 
one has g = 1.31, and for larger values of k, g will be even larger, meaning a 
more reliable predictions for the mean field method. 

The low energy field theory describing a system tuned to be near point P, 
will have a very interesting particle content. It will consist of gapless fermions 
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living on one surface in momentum space, massless fluctuations in the phase 
of the condensate and massive but very hght fluctuations in the magnitude of 
the condensate. It would be interesting to find some observable that might be 
experimentally measured in order to probe such a spectrum. 

The fact that the mass is particularly small also implies that quantum cor- 
rections may significantly alter its value. A rcnormalization group analysis with 
these three degrees of freedom can clarify how beyond mean field corrections 
may shift its value. We leave this for future work. 

But even before such a detailed study, we propose a striking consequence 
of our results. The correlation length tq, or the typical length scale at which 
the magnitude of A varies in field configurations that arise when the system is 
excited, is inversely proportional to the mass of the Higgs mode of the system. 
For example, rg governs the size of the outer core of a vortex configuration in 
a superfluid phase. This can be seen more concretely by writing the classical 
field equations for a condensate of form A(r) = (A + p{r))cxp{iip{(j))), where 
A is the ground state value of the condensate and (r, (/>) are the cylindrical 
polar coordinates with the vortex at r = 0. (See [38[ for reviews and references 
therein.) The boundary conditions for the field p{r) arc that p should tend 
to —A at the center of the core (where the small fluctuation approximation 
begins to break down) and should tend to as r tends to infinity. For a vortex 
configuration, (p winds around by a multiple of 27r as we traverse a loop around 
the vortex. Sufficiently far away from the inner core of the vortex, the spatial 
derivative of f^fj)) does not contribute significantly to the equation of motion, 
and the classical field equation for static p(r) is, 

p(r) - rgvV(r) == const. , (41) 

where 

ro = VE/i3C) (42) 

From Eq. (J4T|) . it is clear that p will decrease from a value close to to a value 
close to —A, as we go closer to the inner core of the vortex, over a length scale 
tq. The fact that C is numerically small close to the point P in parameter 
space (Fig. [3]), will manifest itself in an increased size for the outer core of the 
vortex. 

This is admittedly a simplified discussion. For example, to construct an ac- 
tual vortex solution, it will be important to include the 77^ term in the effective 
action. But the coefficient of this term is dimensionless, and would not intro- 
duce any additional length scale in the problem, and hence we expect our basic 
argument to remain valid in such a detailed study [56[. 

To see the effect quantitatively, we plot in Fig. ([5]), the outer vortex radius 
square, Tq, as a function of a position away from the center of a harmonic 
trap. We use a standard harmonic trap which models a potential in optical 
lattices. The trap parameters used are u = 1.25 x 10~^'^eV, which gives for 
m = 5.61 X lO^eV for Li, a potential muj'^r^/2 = u}{r/rQf/2 with ro = ST.SeV"^ 
. At the center of the trap ^ = -8 x 10^'^eV. The splitting 5p = 1.15 x lO^^eV 
is constant throughout the trap. i?trap, the distance from the center at which 
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Figure 5: (color online) Outer vortex radius as a function of the radial position r in a spher- 
ically symmetric trap at vanishing temperature. Distances are scaled by iJtrap, which is the 
distance at which the effective chemical potential becomes equal to —l/(2ma^), where a is 
the scattering length. The trap parameters are given in the text. On the right of the vertical 
green line the excitations are gapless. On the right of the vertical black line the system is 
in the normal phase. The value of the vortex outer radius changes only very slowly in the 
gapped region and increases monotonically in the gapless region. The divergent large value 
of the outer core vortex radius corresponds to a transition to the normal state where vortex 
does not exist. Notice that the vortex radius exhibits change in the derivative at the gapless 
point which may serve as a signature for the gapless phase. 



at which the effective chemical potential becomes equal to ~l/{2ma'^) is then 
4.6 X lO^eV"^. We scale the distance by this radius. We want to be in the BEC 
side and choose a = 1 x lO^'^cV"^. With these parameters, the gap at the center 
of the trap is 8.29 x lO^^eV. The trap parameters are chosen as an illustration 
of what effects can be seen by choosing a trap which has a substantial volume 
in a gapless phase. 

As we go out from the centre, the effective chemical potential fi — V{r), and 
therefore A decreases and at r/i?trap ~ 0.12 (corresponding to the vertical green 
line in Fig. [5]) we move into the gapless regime. In the gapless region the radius 
of the vortex increases monotonically until it formally diverges as we enter into 
the normal state with no superfluid vortices. 

The increase of the radius of the vortex core with increasing mismatch in 
the gapless region can be qualitatively explained comparing the kinetic energy 
of a superfluid element close to the superfluid vortex with the "condensation 
energy" associated with the superfluid phase; the condensation energy being the 
difference between the free energy in the homogeneous phase and in the normal 
phase. The definition of the vortex radius is by itself ambiguous, because there is 
no abrupt transition from the superfluid phase to the normal phase and various 



definitions have been proposed, see 6.17. [57[. However the length scale at which 



the condensation energy is equal to the kinetic energy should give a qualitatively 
correct result. In particular we expect that the vortex radius estimated with 
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this methods should increase steeply in the gapless phase. The kinetic energy 
of a fluid element close to a vortex is given by 

Ek = n\mv^ = -^ , (43) 

where n is the local superfluid density, m is the mass of the atom and the 
velocity of superfluid matter near a vortex is given by 

v(r) = ^e,, (44) 

with r the radial distance from the center of the vortex and e^ the tangent unit 
vector. As we approach the vortex core the velocity increases and consequently 
the kinetic energy increase. In principle the velocity and the kinetic energy 
diverges for r — > 0, signaling that a certain point, i.e. at a certain value of r, a 
phase transition to the normal phase has to take place. 
The condensation energy is given by 

-Econd = necond , (45) 

where Ccond is given by the difference between the free-energy densities of the 
superfluid phase and of the normal phase. Equating Eq.((331) to Ea.(|^5|) we find 
that the vortex radius is given by 



V Smecond ' 

In the gapped region the energy difference between the superfluid phase and 
the normal phase is not strongly dependent on 5^, thus fo is approximately 
constant. In the gapless phase the condensation energy continuously decreases 
on increasing asymmetry and tends to zero at the boundary between the gapless 
and the normal phase. Thus, fg continuously increases in the gapless phase and 
at the boundary between the gapless phase and the normal phase fp diverges. 

Notice that this deflnition of the radius of the vortex core has to be taken 
with care, because in the normal phase the condensation energy vanishes and 
Eq. (1^5)) seems to suggest that fo diverges. However, in the normal phase there 
is no superfluid motion, consequently there are no superfluid vortices. 

The properties of vortices in the gapless region can not be tested in the 
BCS regime however because there is no stable gapless region in the BCS side. 
Therefore we look at the BEC regime. Qualitative explanation of why the vortex 
size increases with the mismatch, Eq. (P5|) , can not be applied to the BEC region 
with negative chemical potential, /^ < 0. Below wc consider another approach. 

We considered the outer core radius, Eq. (j42p . which is obtained from an 
expansion around the nontrivial vacuum state A ^ 0. To further analyze the 
vortex structure, we obtain the inner core radius which uses an expansion around 
A = state, i.e. Ginsburg-Landau expansion. In [Appendix C| we construct 
the Ginsburg-Landau functional to the fourth order, and derive the equation 
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5n/kT 



Figure 6: (color online) Left panel: Numerical solution of the boundary value problem ODE, 
Eq. 1 148 I I. The center of the vortex is at the origin of the axes coordinate and we report 
the plots of the condensate, full line, and its derivative, dotted red line, as a function of the 
distance from the center of the vortex. The condensate saturates at the boundary value, 
A(r — > oo) = 1.4. Parameters are fi = —1, S/i = 0.4, m = 10, a^ = 0.5, where we are using 
units of T. Right panel: Normalized vortex size as a function of the Fermi momenta mismatch. 
The vortex radius was extracted from the condensate configuration by two methods: based 
on the condensate (as it reaches the value A = 0.7), lower full curve, and the condensate 
derivative (as it reaches A' = 0.03), upper dashed red curve. There is a transition to the 
gapless state around Sfi = 0.6. At this point the slope of the upper curve increases. The slope 
of the lower curve changes at this point too, but the curve is smoother. 



obeyed by ?/(r) in a vortex configuration using the Time Dependent Ginzburg 
Landau equation (TDGL) [52|, 



(a + bvirf - ^V^) ,y(r) = 



(47) 



with the boundary conditions ri{r = 0) = 0, ri{r — > oo) = r]Q. The center of the 
vortex is at r = 0, and the expressions for the coefficients a, 6, c are reported 
in Eq. (jC.5[) . Notice that the TDGL equations for the vortex configuration are 
valid for T ^ Tc where the gap is vanishing small. Introducing ri{r) = e^'^f{C)r]o, 
with C = ry^2mT]o a dimensionless variable, we obtain the TDGL equation for 
the radial part of the condensate configuration 



1 d 






+ ~af- bf 







(48) 



with boundary conditions /(O) — and /(oo) = 1, and where the expression 
for the coefficients a, b, c are reported in Eq. (jC.Sp . At nonzero T, we solve 
numerically the TDGL equation for various values of dfi. In the left panel of 
Fig. ini we report the result of the numerical solution of the condensate and of 
its derivative for Sfi/T = 0.4. In the right panel of Fig. ^ we report the value 
of the vortex radius as a function of 5/i. The two curves correspond to two 
different definitions of the vortex radius, based on a certain value of either the 
condensate or it's derivative. Although calculations are done at nonzero T, one 
may assume that the same trend holds for vanishing temperatures. 
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5. Conclusions 

In conclusion, our results show that as wc move into the gapless regime, 
the outer radius of the vortex increases sharply. This rise may be observed in 
experiments done with cold atomic gases trapped in a magnetic trap. 

If the parameters of the trap, namely, the number of particles of the two 
species, Ni and N2, and the scattering length, a, are tuned such that there is 
a sufficiently wide region in position space where the atomic system is in the 
gapless BEC phase, this dramatic effect can be seen. We leave the precise de- 
termination of parameters a, TVi and N2 for future work, but it will presumably 
require very flat traps to realize this phenomenon in a wide enough region in 
the system to be observed cleanly. 

The properties of vortices in the gapless region have been studied previously 
in 58|, |59[ who have concentrated on the interaction between two vortices in this 



regime. The vortex core structure in imbalanced superfluids has been studied 
in [60| who have focused on the occupation number of particles that determine 
the "visibility" of vortices. The authors of [61[ used a Bogoliubov-de Gennes 
approach to solve for a vortex core state in fermion mixtures with unequal 
masses. They found that the vortex core is mostly occupied by the light mass 
fermions and that the core density of the heavy- mass fermions is highly depleted. 
We believe that their study points towards the gapless phases, however their 
calculations arc more involved. Our study provides motivation to study a new 
observable, namely the size of the core of a vortex, in the gapless phase. 

Finally, we comment about the instability toward the formation of a non 
homogeneous phase. One can see Fig. [3] that in the strong coupling regime 
the coefficient B is always positive. Indeed, along the curves corresponding to 
K = 1 and K ~ 1.71 in Fig.[3]the coefficient B is positive and large. This means 
that there is no instability toward a LOFF-like phase. This is consistent with 
the results of Ref. [23i|, where the LOFF phase was found to be favored in the 
weak coupling regime only. Indeed from Fig. [3] one can see that in the weak 
coupling limit, the curves corresponding to k = or k = —0.5 pass through the 
region where B is negative and this indicates that it is possible to have a non 
homogeneous LOFF phase. But to really check the favorability of a LOFF-like 
phase in this region, in a small 77 calculation, one should expand around the 
solution with A = and not A ^ 0. The reason being that the phase transition 
from (some) LOFF phases to the normal phase is second order and one can 
study how fluctuations drive the system from the homogeneous normal phase 
to a non homogeneous phase. We leave such an analysis for future work. 
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Appendix A. Relations between the Meissner mass, the Debye mass 
and the coefficients of the effective action describing 
the phase fluctuations 

We show the equivalence between the screening masses and the coefficients 
in effective action for the Nambu-Goldstone mode. This equivalence can be 
anticipated from gauge invariance if we gauge the global symmetry associated 
with total number conservation J62l | . The gauged action is invariant under a local 
rotation in the phase of the fermion fields "0^, with /3 = 1,2, accompanied by 
gauge transformation on the four vector (^07 A). (We apologize to the reader 
that we use the same Latin character for the gauge fields as well as for the 
operator defined in Eq. (^0]) . They can be easily distinguished because gauge 
fields always appear with a subscript (Aq) or in bold font (A).) The condensate 
in Eq. (j31) spontaneously breaks this gauge symmetry, and therefore by the 
Anderson-Higgs mechanism, the gauge field components A, acquire a Meissner 
screening mass. The Aq component of the gauge field is instead Debye screened. 
(In the gapless regime, there is an additional contribution to the Debye mass 
from the fermions in the blocking regions, that we do not consider here [47[. 
Adding this contribution to the pairing contribution that we calculate, gives 
the net Debye mass square for the system.) In this Section we explicitly show 
that the screening masses can be related to the coefficients that appear in the 
effective Lagrangian describing the Nambu-Goldstone bosons. 

On gauging the quadratic part of the Nambu-Gorkov action in Eq. ([TB]) we 
obtain 

£ = (-01, -02) 

( idt + ^^^^+gAo + fi + Sfi -A V ^M 

(A.1) 

A gauge transformation is given by ■0/3 -^ ipp^^"'"^^ ^or fermions, leading to 
(V'iV'2) oc A(a;) -> Ae^*"^^) for the condensate, and by A ^ A + ^Va, Aq -> 
Aq + -dta for the gauge field. Since the term with (5/i — > ~6fi is common in 
all the following expressions, we will stop writing it explicitly from Eq. (|A.2p 
to Eq. (|A.17|) and carry it implicitly, and only write it in the final expression 
Eq. (IXT81) . 

First we establish the coefficients in the effective action of Nambu-Goldstone 
mode, in the absence of any external gauge fields. These are space-time depen- 
dent phase rotations of the A-field which we parameterize as A — > Aexp(2iQ;). 
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The quadratic piece of the Nambu-Gorkov fields has the form, 



(V'l>2) 



idt + 2^V2 +^l + S^l -Ae^'" \ . Vi 



-Ac-2*" idt - 2^V2 - ^ + ^M y W' 



"-2 



(^r>2)(0 + ^)( JJ ) (A.2) 



where O and y arc given by Eq. ([T8|) . To the relevant order this phase shifts A 
to A + 7], with fluctuations j] being given by 77 = A(2ia). The effective action 
for a is obtained by integrating out the fermions and can be written directly 
from Eq. ([5T|) by substituting this value of 77. We obtain, 

5GoM.to„c--A (^^j L"(fcH-fc)| Dip)D{p + k) 1' ^^-^^ 

where fcg = ?fc4. 

Next we remove the phase from the condensate by redefining the phases of 
the fermionic fields. This will give rise to non-zero values of the gauge fields 
and we evaluate the screening masses of these gauge fields. Redefining ^p = 
tjjp exp (ia) and acting with a derivative operator, the quadratic part of the 
action can be written as 

(J, %— ,7, p-^ f *^* + 2^^" +^^ + 6f^ -Ae2- \ V^ie- . 

^^^ ' ' ^^' ) (, Ae-2- ^^t - ^V^ - M + J^ j ( ^*e-» ^ 

= (V;i*,^2)(0 + f)(|) (A.4) 

where O is given by Eq. (jf 8p . and y is given by 

\ n r) n, -I- IV'V"+'V"-V-I I IV"; 

V " '^*" "•" 2m "'"2m 

(A.5) 
which includes first and second order terms in a. This Lagrangian is exactly of 
the form given by Eq. ()A.ip with g{AQ, A) = {—dta, — Va) and hence the gauge 
boson masses can be read from the Lagrangian describing the a fields. 

To show that the masses we obtain this way are the same as the coefficients 
obtained by treating a as the Nambu-Goldstone field (Eq. (jA.Sp ) we explicitly 
calculate the second order correction to the action using Eq. (fTT]) . We will 
analyze separately the quadratic term that contain only spatial derivatives of 
a, the term that contain only time derivatives of a and the mixed term. 

Consider the quadratic spatial component, (Va)^. Both linear and quadratic 
terms in V contribute to this part of the effective action. Indeed from the 



(V-iVa+iVa-V) 


(Vaf 











2m 


2m 















fta + 


(V-iVa+iVa-V) 
2m 


■ + 


2m 
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expansion of the action we obtain 



1 



-25ipLi =Tr(6-iy) - ±Tr(0-iy6-iF)|9,„=o 



^ - e A{p) - A{p) 



1 ((2p + k) • k)2 A{p)A{p + k) + ^(p)A(p + fc) + 2A2 



2m i?b) 



} 



2 (2m)2 D{p)D{p + fc) 

(A.6) 

Note that differential operator in Va acts on external in- and out-going legs, 
i.e. produces zk and — ik, while the other V acts inside the loop, producing 
ip and i(p + k). Using the definitions of A{p) and A{p) (Eq. ([20])) the term 
proportional to A{p) — A{p) on the right hand side of Eq. (|A.6p simplifies to 

k.p 

The second term on the right hand side of Eq. (jA.6[) can be simplified by 
noticing that 

C2n -u k") • k _ _ 

2m =e(P + k)-^(p)=A(p)-A(p + fc) + fco^A(p + fc)-A(p)-fco- 

(A.8) 
Then the second term on the rhs of Eq. (|A.6p is given by 

-j ga(fcH-fc)|--(e(p)-e(P + k)) D{p)D{p + k) . 



^^'Y-.^.^.. m/ ^2 2(e(p)-e(P + k)) 



'^a(fcH-fc){-A 



2 



yy f^ ^ ' ^ 'V D{p)D{p + k) 

^^(C(\ CI , .,^^ A{p)A{p + k)-A{p)A{p + k)) ^ /Tx2 
- 2fco(e(p) - C(P + k)) j(p)zj(p + fc) 1 - (y j '^ ' (A-9) 



where 

,. . i:o(.)a(-.){ '^'■" -j;„P/ -" (i(p) - .4(p)) - '^'P'^-^^'; + "" (i(p + .) - A(p + .))} 

- E"(^»(-'=){^rap)-«p-^-«p+^)}-i:{-|S^} 



fc,p ^-^' k,p ^^' 

Here we have used Eq. (|A.8p to rewrite the combination ^(p) — ^(p + k) and to 
extract extra power of D{p) = A{p)A{p) — A^ and D{p + k) in the numerator. 
Combining all the terms together, we get for the spatial component of the 
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effective action 

''^P (A.IO) 

^ 1 t lflr.\ ffr. ^ 1.^^ i(p)i(p + fc) - A{p)A{p + fc) ^ 

Now wc turn to the quadratic term containing only temporal derivatives, i.e. 
(dta)"^. Only terms quadratic in V contribute to this part of the effective action 

(Eq. m), 

^A^V- /,N . ,^( HMp)Mp + k) + A{p)A{p + k)-2A^\ (A-11) 



'^a(fc)«(-fc){ 



y/ ^ ' ' ' 'i 2 D{p)D{p + k) i ■ 

We can simplify this expression, noticing that one can write 

fco = A{p) - A{p + k)- (e(p) - e(p + k)) = A{p) - A{p + k) + (e(p) - e(p + k)) 

(A.12) 
whereupon the right hand side of Eq. (|A.lip can be rewritten as 

J. \ X ^ / , \ / , \ \ K '^ D 



(y) E"(^)"(-^){^' 



D(j>)D{p^k) 
^h(i:(r.\ '■'^ ' ^^^^ Mp)Aip + k)-A{p)A {p + k)x fT^2 

(A.13) 



where, 

D{p + k) D{p) 



Sr (i\ ( ;.N^- ( Aij, + k) + A{p + k) A{p) + A{p) ^ 
s = ^a{k)a{~k)-k,[ ^^-^^ _^j=o. (A.14) 

Then, the term of the effective action containing the temporal derivatives of a 
turns out to be given by 



S^o.^=(^ Ea(fc)a(-fc){-A^ 



,y/ f^ ' ' ' 'I D{p)D{p + k) 

'"'" (A.15) 

+ 4fco(e(p) - e(p + k)) D{p)D{p+k) 1 • 

Finally we consider the mixed component, (Va)(9f a). The only contribution 
to this part of the action comes from the quadratic term in V (Eq. ([T71) V 
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Therefore, 

2 
■^A^V- nx I ,xr, (2p + k) • k i(p)i(p + fc) - A(p)A(p + fc) ■ 



2'5i'Ld--oT'-(0^^^0-^^)|„.i.ed 



> a(fc)a(— fc)-^ fco- 



y ) ^L^ ^ > '- 'I " 2™ D(p)D(p + k) 

k,p 



(A.16) 



We use Eq. (jA.8|) to simplify the mixed term, and we obtain the following 
contribution to the effective action 

c(2) f^\^X^ (i,\ I l.^/ ^uici \ CI , , ^^ Mp)A{p + fc) - A{p)A{p + k) ■ 
'5— = (^j Ea(fc)a(-fc)|--fco(e(p)-^(p+k)) ^(p)^(p + fc) , 

(A.17) 
The sum of the three terms, Eq. (|A.10|) . (|A.15|) and (|A.17|) . gives the effective 
action for the gauge field 



(2) , c(2) , c(2) _ A2r^^'Y-.YM.y , J ^0' " (^(P) " ?(P + k))^ 

latial tcinpc 

((5^ — > — (5/i) 



5^^' =Ctial+5i^ipo.al + 5^Ld--A^(^) E"(^)"(-^){ D{p)D{p+k) 

k.p 



(A.18) 

We recognize that by putting fco = in 5^^) -^ve reproduce the Meissner mass, 
and by putting k = 0, i.e. ^(p) — ^(p + k) = 0, we obtain the Debye mass. 
Comparing Eq. (|A.3|) and Eq. (jA.fSp . we see that coefficients in the effective 
action for the Nambu-Goldstone mode coincide with the corresponding screening 
masses as we set out to show. 

Appendix B. CoefRcients of the low energy Lagrangian 

The sum over the Matsubara frequencies can be done analytically noticing 
that if f{x) is a function with no poles then one has 

^ E /r^ = I tanh[£-/(2T)] f{-£) , (B.l) 



lUJn + £ 2 

n— — oo 



for w„ = (2n + 1)ttT. Upon substituting this result in Eq. ([^^ one obtains the 
usual form of the gap equation at non vanishing temperatures: 

where for infinite volume, the sum over p can be replaced by an integral over 
three-momentum p and where we have defined, 

.9(e) = -(tanh[ ^^ ] + tanh[ — 5^^]) ^ nf{-S^i - e) - nf{-6n + e) , 
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with Uf the Fermi-Dirac distribution function. 

After evaluating the Matsubara p4 sums using Eq. (|B.ip , we get the following 
expressions, 

^i(^) = ^fE-{(5(^i)-.9W)(t— ^ , ' ) 

16 V ^-^ eie L ' \ko + ei — e Kq — ei + e/ 

p 

+ (ff(ei)+5W)(— ^ , ^^ ^ )| 

'Vfco-ei-e fco + ei + e/J 



- - e ^'0 — £1 + e 
'V/cQ-ei-e ko + ei+e/J 

where the sum is an integral over the three-momentum p and we have indicated 
with ei and ^i the quantities e(p + k) and ^(p + k) respectively, with e and ^ 
the quantities e(p) and ^(p) respectively, and where fcg = ik^. 

The expressions of the coefficients A, B, C, D, E and F in Eq. ([5^ are given 
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by 



A 



B 






8y ^^e-^ 
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81/Z^l-9,„2g3 
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1 1 



1 1 Y-r ^ + /^ /^ + A'~ 



:/ 4y2^l 



me"^ 



me^ 



c ^ Ai^ni-M 



D = ii-V ^ 

S 1/ Z^ ^ f 5 



8V 



E = 



8V ^V^V e>m 



p2(e'* + 9C2A2-2A4) 



9CA^ p^ 

A I 



9^2 A2 _ 2A4) 



3eA2 p2(-3eA2 + A'i) 
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,2p2C'A2 



-3 1 
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^ /20A-* 5A2 
3m V e' e' 
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37Ti \ e^ 



7A2 



■ar 
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8A4 
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^ Lgm 
llv i 

4y Z^^gS 



5A2 
3A2 
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2^/lOA'' 7A2 
3^V~^6 ^ 

i4 qA2, 
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2^ /4A4 3 A" 



3m V 
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)]} 



(B.4) 



where g' refers to the differentiation of g(e) with respect to e. 

To evaluate the integrals we use the following relations. For any function 
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/(e,^) we have, 
35/ 
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(27r)3 



9"'f = 



^('5/^-A), , 
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'2T2 



'A^m 



(V-A2)(3/2) V^~d^V^ 



P+ 






p+ 



-A2m25/j / / 



(5/i2-A2)2V p_ 



p- P+ 



-P-- 



2m2V / 1 d //ex .^±(f± 



(Sn^ - A2) V ^'" de' ^ 5 



(f) 
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(5/x2 _ A2) I p_ de W ^ P- ^ P+ do ? 



P+ 



3a,,3 



m S^ 
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(P^ 



p+ 



(B.5) 



where p± = 6'((5/i-A)6l(/i± v/V - A^) ^J 2min ± y/ S n^ - A^), e(p+) = e(p_) = 

Sn and C(p+) = -C(P-) = ^(^M - A)^,5/i2_ A2. In Eq. dHS]), all alge- 
braic terms featuring p± appear with a corresponding product of 6 functions, 
9{5ii — A)9{ii ± \j5jj? — A2), which we have omitted for clarity. Whenever p± 
appear as limits of the integrals, we can simply use the definitions of p± given 
above to obtain the correct answer. 

We analyze the values of the coefficients A, B,C,D,E,F at T = 0. This 
can be done by taking the limit T — >■ in Eq. (jB.4p . We obtain the following 
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expressions: 
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(B.6) 



The equation for E is more complicated and hence we do not give the detailed 
final expression which however can be obtained from the following equation: 
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where, 



C/20A4 5A2 2^^ 2^^10A'^ TA^ 1 



3V £6 £4 g2y 3 V £6 

^i^SA^ SA^N^ ^2/i^4A4 3A2 



3 V e5 £3 J + 3 (, £5 g3 



For completeness we give the expression for F which can be evaluated simi- 
larly, 
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Appendix C. TDGL equation for the vortex core states 

To analyze the vortex core structure we derive the Ginsburg-Landau func- 
tional expanding the action 

d^xM^-iTrlogS-\ (C.l) 

around a state with A = up to the fourth order in 77, obtaining 

where A, A are defined in Eq. ((20)) . Wc rewrite the action as 
T 



^^'^+^^'^ = 5E^(-^>?*(^)^2(fc) + (f7(-fc)r(fc))V4(fc), (C.2) 



fc 
where J4 can be written as 

Mk) = AtjY ^^i ^ + {5^l^-5^i)\ . (C.3) 

We perform the fco Matsubara sum, take po — 0, the limit p ^ and obtain 
the effective action 
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where the coefficients are given by 
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where n± = n(^fc ± (5/i), n is the Fermi distribution, n'{x) = ^i Cfe = ^ ~ l^i 
jjL = {jjLi + yLt2)/2, ^A* = (y^i ~ /^2)/2 and wc take /^i > ^2- The equation of 
motion, 5Seff/Sr]{p) = 0, is given by 



(C.6) 



a + bv{pf+cj-]r,{p)^0, 



and in configuration space 

(a + b7irf^^V')rjir) = 0, (C.7) 

V 2m / 

with the boundary conditions ri{r = 0) = 0, ri{r — > 00) =770- At T 7^ 0, we have 



m 



(2toAo)3/2 1 



47ras Aq 

(2mAo)3/2 1 
A^J 16^ 
(2mAo)3/2 1 



167r2 



cl/2dx 



A3 i67rnyo {x-pr 

[n'{x - p+)+ n'{x - p-)]) 

Q (X — P) / 

(2mAo)3/2 1 -^ 



dx „ f°°x^/^dx, 
-—pr - 2 / [1 - n[;x - p+) - n(x - p_ 



[l-n(a;-p+) -n(a::-p_)] 



)l) 



A3 167r2 



(2mAo)3/2 1 / p xi/2da; 



^0 



IGtt^ VJq (x — p)' 



■[1 - n{x - p+) -n(.T-p_)] 



ci/2dx. 



a;-p 
2 /•°° a;3/2dx 



i'(a:;-p+) + n'(a;-p_)] 



3 Jo (2; - P)^ 
3 7o a; - P 



[1 - n{x - p+) - n{x - p_)] - - / 7 ^Wi^ ~ P+) + "'(2^ " 



n"(x-p+)+n"(x-p_)]) 



3 Jo (a; - P)^ 
(2mAo)3/2 1 



^0 



167r2 



(C.8) 
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where we introduced T = Aq, and the dimensionless variables x = fc^/27nAo, 
p = fi/Ao and v = Sfi/Ao; zeros of the quasiparticle energy are ai p± = p±v. 
Introducing rj{r) = e'^'^ f {C)Vo with dimensionless ( ~ ry'2mrio, we obtain 
the TDGL equation for the vortex core at T 7^ 

cdcvdcy c- 



S(717(C3!r)-^)+5/-&/^ = 0, (C.9) 



with boundary conditions /(O) = and /(oo) ~ 1. The coefficients a, &, c are 
given in the right hand side of Eq. (jC.Sp . Coefficients a = b = c = 1 correspond 
to a supcrfluid ideal Bose gas discussed by Landau |63j , where a vortex filament 
has macroscopic thickness. Here we are able to study both regimes of BCS and 
EEC, BCS-BEC transition, as well as the situation with nonzero mismatch, 
Sp, y^ 0. We solve this second order ODE numerically for different values of dp. 
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